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Abstract 

We investigate the phase diagram of a three-dimensional associating gas (ALG) model. This 
model combines orientational ice-like interactions and "van der Waals" that might be repulsive, 
representing, in this case, a penalty for distortion of hydrogen bonds. These interactions can be 
interpreted as two competing distances making the connection between this model and continuous 
isotropic soft-core potentials. We present Monte Carlo studies of the ALG model showing the 
presence of two liquid phase, two critical points and A density anomaly. 
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I. INTRODUCTION 



Water is an anomalous substance in many respects. Most liquids contract upon cooling. 
This is not the case of water, a liquid where the specific volume at ambient pressure starts 
to increase when cooled below T = 4°C 111. Besides, in a certain range of pressures, water 

n n 

also exhibits an anomalous increase of compressibility and specific heat upon cooling [2|]-[4J]. 
Far less known are its dynamics anomalies: while for most materials diffusivity decreases 
with increasing pressure, liquid water has an opposite behavior in a large region of the phase 
diagram j^-jisl]. 

It was proposed a few years ago that these anomalies are related to a second critical point 
between two liquid phases, a low density liquid (LDL) and a high density liquid (HDL) |l4|. 
This critical point was discovered by computer simulations. This work suggests that the 
critical point is located at the supercooled region beyond the line of homogeneous nucleation 
and thus cannot be experimentally measured. Inspite of this limitation, this hypothesis has 
been supported by indirect experimental results [lf| . 

Water, however, is not an isolated case. There are other examples of tetrahedrally bonded 
molecular liquids, such as phosphorus 0, Q] and amorphous silica [ljj, that also are good 
candidates for having two liquid phases. Moreover, other materials such as liquid metals 



201 ] and graphite 



2l| also exhibit thermodynamic anomalies. Unfortunately a coherent and 



general interpretation of the low density liquid and high density liquid phases is still missing. 

What kind of potential would be appropriated for describing the tetrahedrally bonded 
molecular liquids, capturing the presence of thermodynamic anomalies? Realistic simula- 



tions of water 



22|]- 



241 ] have achieved a good accuracy in describing the thermodynamic 
and dynamic anomalies of water. However, due to the high number of microscopic details 
taken into account in these models, it becomes difficult to discriminate what is essential to 
explain the anomalies. On the other extreme, a number of isotropic models were proposed 
as the simplest framework to understand the physics of the liquid-liquid phase transition 
and liquid state anomalies. From the desire of constructing a simple two-body isotropic 
potential capable of describing the complicated behavior present in water-like molecules, a 
number of models in which single component systems of particles interact via core-softened 
(CS) potentials have been proposed. They possess a repulsive core that exhibits a region of 
softening where the slope changes dramatically. This region can be a shoulder or a ramp 
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25|]-[43|. Unfortunately, these models, even when successful in showing density anomaly and 
two liquid phases, fail in providing the connection between the isotropic effective potential 
and the realistic potential of water. 

It would, therefore, be desirable to have a theoretical framework which retains the sim- 
plicity of the core-softened potentials but accommodates the tetrahedral structure and the 
role played by the hydrogen bonds present in water. A number of lattice models in which the 



tetrahedral structure and the hydrogen bonds are present have been studied 
them, is the three-dimensional model proposed by Roberts and Debenedetti 



441- 



One of 



46j and further 

studied by Pretti and Buzano [47[ defined on the body centered cubic lattice. According 
to their approach, the energy between two bonded molecules rises when a third particle is 
introduced on A site neighbour to the bond. Using a cluster mean-field approximation, they 
were able to find the density anomaly and two liquid phases. 

In order to investigate the mixtures of water with other chemical species as that present 
in a number of biological and industrial processes, it would be interesting to have a simpler 
model capable of capturing the same essential features observed in water and also being 
able to bridge the gap between the realistic models for water and the isotropic softened-core 
potentials. 

Thus, in this paper we investigate a three-dimensional associating lattice-gas model that 
can fulfill both requirements. Our model system is a lattice gas with ice variables 
which allows for a low density ordered structure. Competition between the filling up of 
the lattice and the formation of an open four-bonded orientational structure is naturally 
introduced in terms of the ice bonding variables, and no ad hoc introduction of density or 
bond strength variations is needed. In that sense, our approach bares some resemblance to 



49) 



that of continuous softened-core models 



51 



52( | . Using this simple model we are able 



to find two liquid phases, two critical points and the density anomaly. The remainer of the 
paper goes as follows. In sec. II the model is introduced and the simulation details are 
given. Sec. Ill is devoted to the main results and conclusion ends this session. 

II. THE MODEL 

We consider a body-centered cubic lattice with V sites, where each site can be either 
empty or filled by a water molecule. Associated to each site there are two kinds of variables: 
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FIG. 1: The two possible states of the water molecules in the body centered cubic lattice. A and 
B molecules are forming an hydrogen bond since two of their arms are pointing each other. 

an occupational variables, <7j, and an orientational one, rf . For n« = the i site is empty, 
and rii — 1 represents an occupied site. The orientational state of particle i is denned by 
the configuration of its bonding and non-bonding arms, as illustrated in Fig 1. Four of them 
are the usual ice bonding arms with r* J = 1 distributed in a tetrahedral arrangement, and 
four additional arms are taken as inert or non-bonding (t* j =0). Therefore, each molecule 
can be in one of two possible states A and B are shown illustrated in Fig. 1. A potential 
energy e is associated to any pair of occupied nearest-neighbor (NN) sites, mimicking the 
van der Waals potential. Here, water molecules have four indistinguishable arms that can 
form hydrogen-bondS (HB). An HB is formed when two arms of iViV molecules are pointing 
to each other with = 1. An energy 7 is assigned to each formed HB. 
In resume the total energy of the system is given by: 



(i,3) 

The interaction parameters were chosen to be e > and 7 < 0, which implies in an energetic 
penalty on neighbors that do not form HBs. From this condition results the presence of 




(1) 
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two liquid phases and the density anomaly. 

The ground state of the system can be inferred by simply inspecting the equation [H and 
taking account an external chemical potential p. At zero temperature, the grand potential 
per volume is Q = e + pp, where p is the water density and e = E/V. At very low values 
of the chemical potential, the lattice is empty and the system is in the gas phase. As the 
chemical potential increases, at p, — —2(e + 7), a gas phase with p = and a low density 
liquid (LDL) with p = 1/2 coexist. In this case, each molecule in the LDL phase has four 
occupied NN sites, forming four HBs, and the energy per site is e = £ + 7. As the chemical 
potential increases even further a competition between the chemical potential that favors 
filling up the lattice and the HB penalty that favors molecules with only four iViV sites 
appears. At p = —6e — 27, the LDL phase coexists with a high density liquid (HDL) with 
p — 1. In the HDL, each molecule has eight NN occupied sites, but forms only four HBs. 
The other four non-bonded molecules are repealed, which can be viewed as an effective 
weakening of the hydrogen bonds due to distortions of the electronic orbitals of the bonded 
molecules. The energy per molecule is then e = 4e + 27. 

Our model may be interpreted in terms of some sort of average soft-core potential for large 
hydrogen-bond energies. The low density phase implies an average interparticle distance 
d<LD = Pld 3 — 2 1 / 3 , whereas for the high density phase we have duo = P~hd = 1- The 
corresponding energies per pair of particles are e^ DL = £+7 and e^ 01 = e+7/2 respectively. 
Figure 2 illustrates this effective pair potential for the case of 7/e = —2. The hard core is 
offered by the lattice, since two particles cannot occupy the same site. 

The system pressure P can be calculated from the grand potential since Q = —P. At the 
g&s-LDL coexistence, P = and at LDL-HDL coexistence point, P = 2e. 

The model properties for finite temperatures were obtained through Monte Carlo simu- 
lations in the grand-canonical ensemble (chemical potential and temperature were kept con- 
stant). The total number of molecules is allowed to change in time by means of the Metropo- 
lis algorithm, where in one time unit (1 Monte Carlo step) we test all lattice sites in order to 
insert or exclude one water molecule. The insertion and exclusion transition rates are written 
as w(insertion) = exp(—A<f)) and w(exclusion) = 1 if A<ft > 0, and w(insertion) = 1 and 
w(exclusion) = exp(+A(p) if A(f> < 0. Here A<p = (3 (e mo i ecu i e - p) - In (2), where e mo i ecu i e 
is the energy of the included (or excluded) molecule, and the factor In (2) guarantees the 
detailed balance. 
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FIG. 2: Effective potential vs inter-particle distance for 7/e = —2. The corresponding energies per 
pair of particles are e I p jDL = e + 7 and e^ DL = e + 7/2 for LDL and HDL respectively. 



The simulations were carried out for lattices with linear size L = 10, and the interaction 
parameters were set to 7/e = —2. Runs were of the order of 10 4 Monte Carlo steps. Some 
test runs were done for L=20, showing no relevantchange in the critical temperatures. A 
detailed study of the model properties and the full phase diagrams was undertaken for an 
L=10 lattice. 

In order to obtain the pressure-temperature phase diagram of the model, the pressure was 
calculated from the simulation data. By numerical integration of the Gibbs-Duhem relation, 
SdT — VdP + Ndp = at fixed temperature, we obtain P(p,T), using the condition that 
P = at p = 0. Since the model presents two first-order phase transitions (from gas to 
LDL, and from LDL to HDL phases), the curves p versus p have two discontinuities and 
hysteresis loops. The hystereses were observed when the simulations were started at different 
initial conditions for a given chemical potential around the transition point. 
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III. RESULTS AND CONCLUSIONS 



The model properties for finite temperatures that were obtained through simulations at 
constant temperature and chemical potential go as follows. For sufficiently low values of the 
chemical potential and at low temperatures, all attempts to insert molecules are frustrated, 
and the total density p remains equal to zero. By increasing the value of p, the molecules 
begin to enter in the system, increasing p and leading to two first-order transitions, one 
between the gas and the LDL phases and another between the LDL and HDL phases. The 
dependence between p and the reduced chemical potential p = p/e for some temperatures 
is illustrated in Fig. 3a. Similarly the number of hydrogen bonds per site is illustrated in 
Fig. 3b. The transition between one hydrogen bond per site to two hydrogen bonds per site 
occurs at the LDL-HDL phase transiton. The coexistence of the gas and LDL phases and 
the LDL and the HDL phases were then obtained from this data. 

In order to confirm the loci of the coexistence lines and the critical points, the histograms 
of the densities were collected during a simulation run. The histograms for four different 
temperatures and chemical potentials are shown in Fig. 4 for illustration. Near a first- 
order phase transition, mainly inside the metastable region, the histogram is double-peaked, 
and the system density fluctuates around two characteristic values. One can obtain the 
coexistence lines by finding the chemical potential in which both peaks have the same height 
(Figs. 4a and 4c). As the temperature approaches the critical value, the peaks converge to 
a single one, and a homogeneous phase appears (Figs. 4b and 4d). 

In Fig. 5 we exhibit the reduced pressure, P — P/e, versus density isotherms. The gas- 
LDL and LDL-HDL first-order phase transitions are evidenced by the presence of plateaus 
in the P.vs.p curves at low reduced temperatures, T = T/ks- 

The plot of density versus reduced temperature at constant pressures shows that an 
inversion of the behavior of density as a function of temperature takes place at intermediate 
pressures, in the LDL phase. At smaller pressures, density decreases with temperature, 
whereas at higher pressures, density increases with temperature. This yields a temperature 
of maximum density for a fixed pressure, TMD, in the higher range of pressures, which we 
illustrate in Fig. 6. 

The pressure-temperature phase diagram is illustrated in Fig. 7. The gas, LDL and HDL 
phases are shown together with the two coexistence lines, the two critical points and the 
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FIG. 3: (a) Density isotherms vs. reduced chemical potential for different temperatures, (b) 
Number of bonds per site vs. reduced chemical potential for different temperatures, p is given in 
units of lattice space and the temperature is in units of ks- 

line of temperature of maximum density (TMD) as a function of pressure. Density versus 
reduced temperature illustrating the two coexistence regions and the two critical points are 
shown in Fig. 8. 

We have shown that it is possible to incorporate some of the microscopic properties of true 
water molecules into a very simple minimal model that still contains some of the ingredients 
of real water without having its whole complexity. 

The model includes orientational and occupational variables, and guarantees the local 
distribution of hydrogens on molecular bonds, without the need of increasing the volume 
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FIG. 4: Histograms of the total density p. (a) the coexistence between LDL and HDL phases at 
T = 0.8 and (b) an homogeneous phase at T = 1.0 near the LDL — HDL critical temperature, 
(c) the coexistence between gas and LDL phases at T = 1.2 and (d) an homogeneous phase at 
T = 1.4, near the gas — LDL critical temperature. 



artificially or introducing artificial orientational variables 45] . In spite of the absence of 
an orientational order-disorder transition [53], the model presents liquid-liquid coexistence, 
with slightly positive slope in the pressure-temperature plane, accompanied by a line of 
maximum density on the low density side, a feature expected for real water. Besides, this 
study points out to the fact that the presence of a density anomaly, with the thermal 
expansion coefficient a < 0, on the low temperature side, and as a consequence, (§^)t > 0, 
does not imply a negative slope of the liquid-liquid line, contrasting with the results for 
most studies of metastable liquid-liquid coexistence in models for water, which suggest a 



transition line with negative gradient 



26) 



The presence of both a density anomaly and two liquid phases in our model begs the 
question of which features of this potential are responsible for such behaviour. Averaged 
over orientational degrees of freedom, our model can be seen as some kind of shoulder 
potential, with the liquid-liquid coexistence line being present only for a repulsive "van der 
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FIG. 5: Reduced pressure as a function of the total density p for some values of T. 
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FIG. 6: Total density as a function of reduced temperature at constant values of the reduced 
pressure. The maximum in the curves give the temperature of maximum density for a given 
pressure. 
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FIG. 7: Reduced pressure versus reduced temperature phase diagram. The gas-LDL and LDL- 
HDL coexistence lines, the two critical points and the TMD line are shown. 
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FIG. 8: Density as versus reduced temperature illustrating the two coexistence regions, the two 
critical points and the TMD line. 
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Waals" potential. The same was indeed observed for continuous step pair potentials [30I. 
for which, however, the density anomaly is absent. On the other hand, a density anomaly 



has been observed in a number of shoulder- 
is the competition between two scales 



31 



ike 



32 



attice models in which the major ingredient 



331 ] . This feature is present in our case. 



Therefore it seems that the competition between two scales is the major ingredient that 
warrantlES the presence of the density anomaly. If, in addition, the model has an attractive 
interaction, two liquid phases and two critical points emerge. 
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